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A directed avalanche model with a control parameter is introduced to describe the transition 
between cohesive and noncohesive granular material. The underlying dynamics of the process can be 
mapped to interface growth model. In that representation, a continuous phase transition separates 
the rough phase and the flat phase. In the avalanche formulation, this corresponds to a transition 
from deep to shallow avalanches. The scaling exponents of the avalanches indeed follow those of 
the underlying interface growth in both phases and at the transition point. However, the mass 
hyperscaling relation is broken at the transition point due to the fractal nature of the avalanche and 
a hierarchy of critical directed percolation processes. 
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I. INTRODUCTION 

Granular avalanches have received much attention 
since sandpile models are used as paradigms of so-called 
self- organized criticality Q]. However, observations of 
critical-type distributions of avalanches in real physical 
systems are still rare, with as a notable exception the 
recent rice pile experiments by Frette et al. It was 
suggested by Christensen et al. || that the anisotropy 
in the rice grains allows more stable packing configura- 
tions in a granular pile, and that this could be responsible 
for the successful observation of criticality. Some of the 
recent attention has been drawn to avalanches in cohe- 
sive granular materials with the premise that cohesion 
will also allow the sand more packing configurations and 
thus increase the likelihood of observing critical scaling 
behavior. While the goal of finding criticality in cohesive 
sandpiles remains to be fulfilled even after the experimen- 
tal work by Quintanilla et al. the effect of cohesion 
in granular avalanches represents an interesting direction 
for a theoretical study. 

In this article, we'll use the discrete-height version of 
the sandbox (DHSB) model introduced in Ref. || for an 
unloading sandbox (Fig. 0) to understand the effects of 
cohesion in directed avalanche systems. In the following 



section, we'll discuss how we can model cohesiveness in 
avalanche systems. In Sec. [II, we'll review the DHSB 
model and introduce a cohesion parameter. Previous re- 
sults in Refs. || || represent a special case of the model 
where the system is in the deep avalanche phase with 
the cohesion parameter p — 1/2. In Sec. [TV], we de- 
scribe the step-flow random-deposition (SFRD) interface 
growth model which underlies the DHSB model and the 
directed percolation (DP) roughening transition of the 
SFRD model. In Sec. [vj we focus on the two determin- 
istic limits of the model and present the exact solution 
in one of these limits. In Sec. VI, numerical results for 



the avalanches in the flat phase of the interface model 
are presented. In Sec. VII, we investigate the scaling be- 



havior at the transition point where the interface rough- 
ness increases logarithmically in time. We show that 
the avalanche-scarred sand surface, while being rougher 
than nonscarred ones, retains the same scaling exponent 
of the roughness in the thermodynamic limit. However, 
we'll also show that at the transition point, the viola- 
tion of mass hyperscaling relation spoils the reduction to 
two independent exponents es tabli shed in Ref. ||. We'll 
summarize our results in Sec. VIII. 



II. TUNABLE PARAMETER FOR COHESION 
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FIG. 1: A sandbox system. The rectangular box is filled 
with sand. One of the retaining wall can be lowered slowly to 
let out the sand in a sporadic way forming distinct avalanche 
events. 



One interesting characters of cohesion in sand is that it 
possesses hysteresis behavior. Consider building a sand 
castle on a beach. It's common sense that we'll need to 
add water to the sand before we can shape it into a stand- 
ing castle. However, without disturbance, the sand castle 
can somehow maintain its shape even after it dries out 
. The moisture in sand increases the cohesion between 
the sand particles J§| and allows one to manipulate the 
sand into a stable shape that, while not as attainable, is 
more or less an equally valid stable shape for dry sand. 

In accounting for this standing-sand-castle effect, we'll 
use the same stability condition for all cohesiveness of the 
sandbox. While, in reality, the space of possible stable 
configurations for wet and dry sand should not be exactly 
identical, in this article, we shall ignore this distinction 
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FIG. 2: The lattice structure of the two-dimensional discrete- 
height sandbox model. It corresponds to a top view of the 
sandbox with the lowering wall located at the bottom. 

to avoid complicating the rules too much. 

On the other hand, the way an unstable sand surface 
topples surely depends on the cohesiveness. In the DHSB 
model discussed below, there are only two possible final 
stable states for any toppling site. We'll call them the 
minimal stable state and the maximal stable state. These 
two states are similar to the angle of repose and maxi- 
mal stable angle in a real sandpile. However, in sandbox 
model, these states are microscopic while the "angles" of 
a real sandpile are macroscopic. We'll use a parameter 
p, which is a real number between and 1, to represent 
the strength of cohesion. In the model, p is the proba- 
bility for a toppling site of the sandpile to settle into the 
maximal stable state instead of the minimal one. For wet 
sand, the p is large, and for dry sand, the p is small. 

III. DISCRETE-HEIGHT SANDBOX MODEL 

With the discussion of the previous section in mind, 
let's review the dynamic rules of the discrete-height sand- 
box model. The surface of a sandbox (see Fig. 0) is 
represented by an integer height variable h defined on 
a two-dimensional square lattice which is tilted at 45° 
with respect to the lowering wall as illustrated in Fig. |^. 
This is equivalent to considering only the lattice points 
whose integer x and y coordinates satisfy the condition 
that x + y is an even number. The lowering wall that 
drives the system by creating unstable sites is located at 
the y = row and the activities in the system propagate 
only in the positive y direction. In our numerical simula- 
tions, the system is periodic in the x direction, which is 
parallel to the driving wall. The sizes of the system in the 
x and y directions are denoted by the numbers of sites 
L x in each row and the number of rows L y respectively. 

As in most sandpile processes, the dynamics of the 
sandbox model is defined by a stability condition, a top- 
pling rule, and a driving method. They are as follows. 
The stability condition of the DHSB is given by 

h(x, y) < min [h(x — 1, y — 1), h(x + 1, y — 1)] + s c (1) 
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FIG. 3: A typical configuration of the discrete-height sand- 
box model before (left) and after (right) a system spanning 
avalanche. Sites participated in the avalanche are shaded 
darker. The system size L x x L y is 32 x 64. 

slope. The unstable sites in the system topple with the 
rule 

h(x, y) — » min [h(x — 1, y — 1), h(x + 1, y — 1)] + ?/, (2) 

where r\ = with probability 1 — p and rj — 1 with 
probability p. (In the earlier studies || [|, the value of p 
is always 1/2 .) This is the only place in the dynamics of 
the DHSB that the cohesion parameter p comes into play. 
The lowing wall which drives the system is implemented 
in the model by randomly picking one of the highest sites 
(xi, 0) on the y — row and by reducing its height by 1: 

h(xi,0)->h(x it 0)-l, (3) 

where i is the Monte Carlo time, which also serves as an 
age index for the avalanches. 

A typical configuration of the DHSB before and after 
an avalanche is shown in Fig. [?| Since the toppling of a 
site on a given row y only affects the stability of the two 
sites immediately above it at the y + 1 row, we choose 
to update the system in a row-by-row fashion. For each 
avalanche, the entire system is stabilized by such a single 
sweep of topplings from y = to y = L y . 

IV. UNDERLYING INTERFACE DYNAMICS 

The underlying interface dynamics of the sandbox 
models is given by the step-flow random-deposition 
(SFRD) models with a two-step growth rule || |6| as 
illustrated in Fig. || The mapping between the sandbox 
system and the interface growth model involves identify- 
ing the y coordinate of the sandbox model with the time 
t of the interface growth. Each stable sandbox surface 
thus can be viewed as a space-time world-sheet config- 
uration of the interface growth. Models similar to this 
generally belong to the Kardar-Parisi-Zhang (KPZ) uni- 
versality class H with the critical exponents a = 1/2, 
(3 = 1/3, and z = a/ (3 = 3/2 which characterize the 
scaling of interface roughness 



with s c = 1, which represents the local maximal stable 



W 2 = (h-h) 2 . 



(4) 
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FIG. 4: Two-step growth of the discrete-height step-flow 
random-deposition interface growth model; (a) Steps flow by 
one unit to the right (left) when its size Aft is negative (pos- 
itive); (b) Each site increases by one unit with a probability 
P- 



Starting from a flat interface at t = 0, the interface grows 
rougher with 



W ~ f 



(5) 



And, after a characteristic time t c ~ L z , the roughness 
will saturate with a value 



W ~L Q 



(G) 



depending on the system size L. 

From the mapping introduced in Ref. Q , the avalanche 
exponents are given by 



u — 1 — Oi 

n = = 2, 



(7) 





FIG. 5: Scar (edge lines of avalanche clusters) configurations 
of DHSB avalanches at the two deterministic limits; (a) p — 0; 
(b) p = I . 



However, the discrete-height version of the SFRD 
model undergoes a DP roughening transition at p = 
p c w 0.294515 similar to those studied by Kertesz and 
Wolf @ also Alon et al. fyj. The KPZ scaling behavior 
only applies when the value of the control parameter p is 
greater than the critical value p c . Below this transition 
point the interface is in a trivial flat state, where, for a 
stationary interface (interface time y — > oo), the density 
of sites at the bottom h = ho layer is finite. The interface 
is thus pinned at this level and its growth rate becomes 
zero. 

At the transition point p — p Cl we find the roughness of 
the SFRD interface diverges only logarithmically in time 



r w =a-z-a=-, 



and 



(J — 1 — z 



T6 



= 4 



(8) 



(9) 



W 2 



(lnf) 7 , 



(11) 



with the exponent 7 ~ 1 similar to that of the Kertesz 
and Wolf's model as well as the restricted version of the 
models by Alon et al.. 



for the distribution functions, Pi (I) ~ l Tl , P w (w) ~ w Tm , 
and Ps (S) ~ 5 TS , of avalanche length /, width w, and 
depth 6. As defined in Ref. §, the avalanche length Z 
(width it; ) represents maximum y (x) distance of the top- 
pling sites from the triggering point while the avalanche 
depth S is the maximum height change of the toppling 
sites. The a in these expressions was eliminated with the 
mass hyperscaling relation 



2a. 



(10) 



obtained from the compactness of the avalanche clusters, 
i.e., assuming m ~ IwS. 



V. DETERMINISTIC LIMITS 

In the two limits, p = 1 and p = 0, the toppling 
process of the avalanches becomes deterministic and the 
sand surface topples down layer by layer. The only ran- 
domness in the process comes from the driving method 
(||), i.e., that we randomly lower one of the highest sites 
at the y = row to trigger an avalanche. The typi- 
cal avalanche scar configurations at these two limits are 
shown in Fig. ^|. These are the edges of avalanche clusters 
left on the surface, some of which are partially erased by 
newer avalanches. 
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FIG. 6: The domains of odd (shaded region) and even (light 
region) sites on a DHSB surface at p = 0, separating them 
are domain walls that no avalanche will penetrate at this de- 
terministic limit. 



A. Domain walls at p = 

The p = limit runs into the complication that in the 
bulk of the system (y > 0) the sand surface goes down by 
2 units at a time. Since Ah = h(x, y) — n\m\h(x — l,y — 
1), h(x + l,y — 1)] = 1 is stable according to the stabil- 
ity condition (|l|), and the sites on the y = row always 
goes down by 1 unit each time according to the driving 
method (||), the sites on the y = 1 row will only topple 
when their heights are 2 units higher than the trigger- 
ing sites and they alway go down by 2 units to the same 
height of the triggering site according to the toppling rule 
(||). All the sites at higher rows will be locked into the 
same even-oddness as the sites triggering their toppling. 
Therefore, after all sites have participated in at least one 
avalanche, their even-oddness will be fixed for all subse- 
quent topplings. This means the even-oddness of a site 
is preserved by the toppling process, and that the lines 
separating the even and odd sites thus form impenetra- 
ble domain walls for the avalanches (see Fig. ||). This 
hinders the applicability of the same type of analysis as 
presented below for t he p = 1 limit. However, the nu- 
merical results in Sec. |v| will show that the same scaling 
exponents as those of p = 1 case control this limit, too. 



B. Exact solution at p = 1 

The p = 1 limit has a nice solution. Since the sites in 
the bulk topple from Ah ~ 2 to Ah — 1, the sand surface 
indeed goes down only one layer at a time without the 
complications as the p = case. An exact solution can 
be obtained by considering the avalanches taking place 
in such one single layer. For a brand-new layer, the two 



boundaries of the first avalanche open up linearly until 
the avalanche spans the system in the x direction and 
leaves two scar lines on the surface. The two boundaries 
of the second avalanche expand until they meet the scar 
lines created by the first avalanche. Then, they turn and 
follow those scar lines until they meet with each other and 
terminate the avalanche. Subsequent avalanches follow 
the same scenario. The maximum distance an avalanche 
cluster can expand from its triggering point to each side 
in the x direction is exactly half the distance from the 
nearest triggering point of the previous avalanches in the 
same layer on that side. As the triggering points are 
chosen in an uncorrelated manner, the maximum width 
w of an avalanche should follow the Poisson distribution 



«»! 



(12) 



if /i is the average distance between the triggering points 
of the previous avalanches in the same layer in the sta- 
tionary state. The avalanche under consideration could 
be any one of the avalanches happening in the same layer. 
Thus, we need to average over the number of avalanches 
n taking place before this one in the same layer. For a 
system of transverse size L x , n — L x /fj,, the integral can 
be carried out explicitly and gives 



f 

Jo 



H w e~v ,1 (w -2)! 2 

; d- = - r^- ~ IU , 



which results in 



n = t w = 2. 



(13) 



(14) 



The same results can also be derived from Eq. (J7]) and 
(0) by assuming z — 1 and a = 0. Since the avalanches 
are compact, the hyperscaling relation ( fl0|) and other 
exponent relations from Ref. ||] hold. 



VI. SHALLOW-AVALANCHE PHASE 

Below the transition point, the underlying interface 
model is in a flat phase where the bottom layer perco- 
lates with finite density. All the information of the initial 
configuration of the interface (the y = row next to the 
wall) is wiped out at a time scale proportional to the 
sizes of the islands higher than the bottom layer in the 
initial state. (Without deposition, the sizes of these is- 
lands decrease linearly in time.) While the underlying 
interface model is in a trivial phase, much like the uncor- 
related stationary state in Dhar and Ramaswamy's di- 
rected sandpile model |L2| , the avalanche distributions of 
the system may still exhibit power-law scaling. The nu- 
merical values of the scaling exponents shown in Fig. |?] 
confirm the power-law scaling of the distributions and 
they are similar to those values found at the p — 1 fixed 
point following z = 1 and a = 0. While an exact so- 
lution is not available in this phase, we can understand 
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FIG. 7: Finite-size scaling (FSS) estimates of the scaling ex- 
ponents versus inverse width (1/w) of avalanche clusters for 
the DHSB avalanches in the shallow-avalanche phase (mea- 
sured at p = 0.1). They are consistent with a = and z = 1. 

the scaling exponent z = 1 from the perspective that the 
DP clusters triggered from single seeds in the percolating 
phase open up linearly I ~ w; and also that roughness ex- 
ponent a — comes from that the interface is flat. How- 
ever, a difference is that while p < p c represents an entire 
phase of shallow avalanches which should be controlled 
by an attractive fixed point, the p = 1 fixed point is un- 
stable in the sense that the scaling behavior falls back to 
the KPZ universality class for any small deficiency in the 
cohesiveness p from the value 1 . 

VII. DP ROUGHENING TRANSITION 

At the transition point p = p c , the interface roughness 
diverges logarithmically thus the (3 and a exponents, de- 
fined by Eqs. (|^) and (g), are both zero. Nonetheless, the 
dynamic exponent z has a nontrivial value zdp ~ 1.582 
originating from the DP nature of the bottom-layer dy- 
namics. Moreover, at the transition point, the avalanche 
clusters lose their compact shapes (see Fig. ||) and we 
should not expect the exponent relations (f7[) — (JlC|) , nor 
the calculation in Ref. for the corrections to scaling 
to remain valid. In this section we will demonstrate the 
break down of mass hyperscaling relation ( ^p| ) and how 
the avalanches affect the roughness of the sand surface. 

A. Breakdown of mass hyperscaling 

At the transition point, the bottom layer of an 
avalanche cluster follows the critical DP dynamics. 
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FIG. 8: Typical large avalanche cluster for DHSB (a) at the 
DP transition point; (b) in the deep avalanche phase (p — 
0.5), triggered at lowering wall boundary at the bottom. The 
length / and width w of each avalanche are as labelled. Black 
area in the cluster of (a) is of sites that topple to the lowest 
height ho of the bottom layer. It shows the percolation of 
the bottom layer. One sees that the avalanche maintains a 
compact structure in the deep avalanche phase while becomes 
more fractal-like at the transition point. 

Therefore, we should expect from the fractal DP clus- 
ter shape that the density of sites at the lowest h = ha 
level goes to zero in the thermodynamic limit for large 
avalanches. However, the overall shape of an avalanche 
consists, in addition, of sites at ho + 1, ho + 2, . . . levels. 
The higher-level sites that participate in the avalanche fill 
into the holes and voids next to the bottom layer cluster 
and more or less bring the avalanche cluster back to a 
compact shape. We can verify this compactness of the 
avalanche cluster by a direct measurement of the ratio 
a/(lw), with a being the area of (or, the number of sites 
participating in) an avalanche. The result is shown as 
the solid line in Fig. ||. The approach to a finite value 
on the vertical axis demonstrates the compactness of the 
avalanche clusters by the existence of a finite area density 
sa 0.2 in the thermodynamic limit. The FSS estimates are 
plotted against 1/ In y instead of 1 jy since the roughness 
of the surface diverges only logarithmically in y, which 
will be elaborated later. 

Contrary to a finite area density, as also shown in 
Fig. [)[ the mass density m/{lwS) (the dashed line) goes 
to zero in the thermodynamic limit. The absence of a 
finite mass density breaks the scaling 

m ~ Iwd, (15) 

which lead, in Ref. ||, to the mass hyperscaling relation 
( |To| ) . The plot of the combined exponent a — z — 2a in 
Fig. |Io| shows the violation of Eq. jl0| ) as the FSS esti- 
mates approach w 1.72 which is much lower than the ex- 
pected value 2 for compact avalanches obeying Eq. Jili|). 
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FIG. 9: The FSS plot of the area density a/(lw) (solid line) 
and the mass density m/(lw8) (dashed line) versus inverse the 
length logarithm (1/lni) for the avalanche clusters at the DP 
transition point. While the area density converges to a finite 
value at the thermodynamic limit, the mass density converges 
to 0. 
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FIG. 10: FSS estimates of the scaling exponents derived 
from the avalanche exponents n , t w , ts for the discrete-height 
sandbox model versus the inverse width (1/w) at the DP tran- 
sition point. The z exponent is consistent with dynamic expo- 
nent of DP universality class zdp — 1.582. The combination 
a — z — 2a < 2 indicates a violation of mass hyperscaling 
relation (JToj) . 



Also shown in Fig. |l0| are the plots for the a, z, and a 
exponents. They are consistent with z — zdp and more 
or less with a = 0. This confirms that the scaling behav- 
ior of the avalanches follows those of the SFRD interface. 
The slow convergence of a is to be expected from the 
logarithmic divergence of the interface roughness. 
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FIG. 11: (a) The roughness of a stationary DHSB sur- 
face(dotted line) compared with the roughness of the SFRD 
model (solid line) versus the double logarithm of time t at the 
DP transition point. The iterated avalanche process makes 
the surface rougher. The dashed line shows the difference 
AW 2 = Wdhsb — Wsfrd between the roughness of the two. 
(b) FSS of the 7 exponents of the logarithmic scaling for the 
SFRD roughness Wi Ft{D (solid line) and the difference AW 2 
(dashed line), both assumed to have the scaling form (Int) 1 , 
versus the inverse of the logarithm of time. In the t —* 00 
limit, AI4 72 scales with a smaller 7 exponent than that of 
W 2 . 



B. Interface roughness 

The remaining question is how the scaling behavior of 
the roughness is changed by the iterated avalanche pro- 
cess. We approach this by looking at the change of the 
global surface roughness itself and by comparing the scal- 
ing of this change to the scaling of the original interface 
roughness. The same analysis was performed in Ref. Q 
which concerns only the p — 1/2 case of the DHSB, and 
it was found that the change in the global roughness by 
the avalanche process only represents large corrections 
to the KPZ scaling behavior of the surface. However, at 
the DP transition point, the interface roughness diverges 
only logarithmically. This makes the scaling of interface 
roughness more likely to be overwhelmed by the change 
in the roughness due to the avalanche process, and we 
generally would not expect the values of the scaling ex- 
ponents to remain the same. In the following, we'll show 
the scaling of the interface does follow the same logarith- 
mic divergence. 

We perform a direct measurement of the global inter- 
face roughness at the transition point. The results are 
shown in Fig. [ll](a). As in the p — 1/2 case, the surface 
is made rougher by the iterated avalanches. The increase 
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in the roughness AW 2 scale as (lni) 7A with the exponent 
7a ~ 0.4 which is shown as the dashed line in Fig. |Tl"|(b). 
Since the interface roughness itself scales as W 2 ~ (lni) 7 
with 7 « 1 which is shown as the solid line in Fig. |ll|(b) , 
the change in the roughness is irrelevant comparing to 
the interface scaling. We can thus conclude that in the 
thermodynamics limit, the stationary surfaces of DHSB 
have the same 7 exponent as the SFRD interfaces. Just 
as the in the deep phase (p = 1/2) of the the avalanche, 
the iterated avalanche process only gives rise to sizable 
corrections to the interface scaling behavior. 



VIII. SUMMARY 

In this article, we introduced the DHSB as a model for 
avalanches in granular materials with variable cohesive- 
ness. This model exhibits a deepening transition from 
a shallow- avalanche phase where avalanches only involve 
a couple of surface layers of the granular material, into 
a deep-avalanche phase where the depths of avalanches 
increase as power laws in their lengths or widths. In 
the deep-avalanche phase, the scaling behavior of the 
avalanches belongs to the KPZ universality class: The 
avalanche clusters scale anisotropically with I ~ w 3 / 2 
and depth increase as S ~ w 1 / 2 . In the flat phase, the 
avalanche clusters scale isotropically I ~ w with finite 
depths. 

In both phases, the mass hyperscaling relation (|To| ) 
based on compactness dlq) of the avalanches holds. On 



the other hand, at the transition point, the hierarchical 
DP structure, pointed out by Tauber et al. Q, for each 
height level breaks this scaling in a subtle way. While the 
mass density m/(lwS) of the avalanche clusters goes to 
zero in the thermodynamic limit, the area density a/(lw) 
remains finite. However, the exact scaling behavior of the 
systems at this DP roughening transition point remains 
unclear even without the iterated avalanche in the DHSB 
model @|l|,|l6|. 

While we are not aware of any experimental study on 
how the avalanche behavior of a system will vary with 
a gradual change in the cohesiveness of the grains, the 
cohesiveness in granular system is known to vary with 
moisture || and grain sizes [l?]]. We thus expect ex- 
perimental studies in this direction to be feasible. The 
DHSB model represents a system with a layered struc- 
ture where the heights are discrete, and the DP nature of 
the deepening transition relies heavily on a well-defined 
bottom layer or minimal stable configuration of the sys- 
tem. It thus wouldn't be a surprise if exact DP scaling 
were not to be observed in the avalanches of most ex- 
perimental sandpiles. Nonetheless, the breakdown of the 
mass hyperscaling relation (10) comes from the fractal 
aspect of the hierarchical DP clusters and is a more fun- 
damental property. It would serve as a hallmark of such 
a transition if it's to be observed experimentally. 
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